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Primordial non-Gaussianity (NG) affects the large scale structure (LSS) of the universe by leaving an imprint on 



the distribution of matter at late times. Much attention has been focused on using the distribution of collapsed 
objects (i.e. dark matter halos and the galaxies and galaxy clusters that reside in them) to probe primordial 
5^ . NG. An equally interesting and complementary probe however is the abundance of extended underdense re- 
gions or voids in the LSS. The calculation of the abundance of voids using the excursion set formalism in the 
presence of primordial NG is subject to the same technical issues as the one for halos, which were discussed 
e.g. in Ref. [51]. However, unlike the excursion set problem for halos which involved random walks in the 
presence of one barrier <5 C , the void excursion set problem involves two barriers S v and S c . This leads to a new 
complication introduced by what is called the "void-in-cloud" effect discussed in the literature, which is unique 
to the case of voids. We explore a path integral approach which allows us to carefully account for all these 
issues, leading to a rigorous derivation of the effects of primordial NG on void abundances. The void-in-cloud 
issue in particular makes the calculation conceptually rather different from the one for halos. However, we 
t-h ■ show that its final effect can be described by a simple yet accurate approximation. Our final void abundance 
\ function is valid on larger scales than the expressions of other authors, while being broadly in agreement with 
■ those expressions on smaller scales. 



H ; 1 Introduction 

A striking feature of the large scale structure (LSS) of the universe which has emerged from studies over the 
last few decades is the presence of a filamentary network or cosmic web in the matter distribution, with galaxies 
distributed along filaments which surround large apparently empty regions termed voids [1, 2, 3, 4, 5, 6]. A 
considerable amount of analytical and numerical effort has gone into understanding the nature of the cosmic 
web (see e.g. Refs. [7, 8, 9, 10], for a review see Ref. [11]). Apart from being intrinsically interesting, the LSS 
has the potential to be a powerful probe of cosmology since it is sensitive to both the expansion history as well 
as the initial conditions (in particular the physics of inflation) . The advent of large galaxy surveys has realised 
this potential considerably over the last decade or so [12, 13, 14], with ongoing and upcoming surveys set to 
significantly increase the precision of the LSS as a cosmological tool (see e.g. Refs. [15, 16, 17]). 

Our focus in this work is on probing the physics of the early universe. In the current understanding of 
structure formation, the present LSS has its seeds in the statistics of the tiny primordial curvature inhomo- 
geneities generated in the very early universe, during the rapidly expanding inflationary phase. The simplest 
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models of inflation involving a single "slowly rolling" scalar field predict that the statistics of these primor- 
dial inhomogeneities are almost Gaussian [18, 19]. Constraining the amount of non- Gaussianity (NG) in the 
primordial distribution therefore provides a unique window into the physics of inflation. Until recently, the 
cosmic microwave background (CMB) radiation was considered the standard tool for constraining primordial 
NG, since inhomogeneities at the CMB epoch are small and the physics can be described by a perturbative 
treatment. In terms of the standard parametrization of the NG, the CMB constraints for the local model 
translate to —10 < f§ c L < +74 [20]. (For a brief introduction to primordial NG see Appendix A.) Over the 
last several years, however, the LSS has been shown to be an equally constraining and moreover complementary 
probe of primordial NG. For example, from Ref. [21] one finds —29 < < +69, already comparable with 
the CMB constraints, with precisions of order A/^£ ~ 10 [22] and A/j^£ ~ 1 [23] being claimed for future 
surveys. These constraints and forecasts rely on the statistics of massive collapsed dark matter halos and the 
galaxies and galaxy clusters that reside in them, chiefly using three tools - the scale dependent bias in the 
galaxy power spectrum [24, 25] , the galaxy bispectrum [26] and the number density of collapsed objects (mass 
function) [27, 28]. For recent reviews, see Refs. [29, 30]. 

On the other hand, the voids found in the LSS are also potentially interesting probes of NG. While any 
given galaxy survey is expected to contain fewer voids than halos, the effect of NG on void abundances is 
opposite to that on halo abundances (see below). Whereas a positive /nl will enhance the abundance of halos, 
it will reduce the abundance of voids and vice versa, making voids a complementary probe of NG. Although 
the literature contains several definitions of what a void exactly means (see e.g. Ref. [31] for a recent review 
and thorough comparison), the basic physical picture is that of a large expanding region which is underdense 
compared to the background. For a review of the structure and dynamics of cosmic voids, see Ref. [32]. From 
the point of view of using voids as statistical probes, detailed numerical studies [33, 34, 35] and analyses [36, 37] 
seem to indicate that the excursion set formalism [38, 39, 40], combined with simplified analytical models which 
track the underdensity of dark matter in a region, provides a good starting point. More precisely, the spherical 
ansatz [41, 42, 43, 36] which we discuss later in a bit more detail has been found to be a very useful tool when 
studying voids. (See also Refs. [44, 45] for a different approach and application.) 

As one might expect from experience with halo abundances, analytical treatments of void abundances are 
subject to some caveats. Firstly, by choosing to describe voids using underdensitics in the dark matter distri- 
bution, one is a step removed from actual observations which involve galaxies. As pointed out by Furlanctto 
& Piran [46], the visually striking empty patches of galaxy surveys are all galaxy voids, and it requires some 
nontrivial analysis to relate their distribution to that of the underlying dark matter. There is also some indica- 
tion from A-body simulations (with Gaussian initial conditions) [34] that analyses of the distribution of voids 
based on the excursion set formalism, such as the one by Sheth & van de Weygaert [37], could fail to capture 
some non-universal features of the void abundance. Another concern is that the spherical assumption is of 
course an idealization - real structure is more complicated. Nevertheless, from the point of view of obtaining a 
better understanding of the underlying physics from an analytical perspective, an excursion set analysis based 
on the spherical ansatz remains the most robust starting point currently available. Once the calculation in 
this simplified setting is under control, one can explore improvements by relaxing the assumptions involved. 
We will adopt this point of view in this work, our goal being to obtain an expression for the abundance of 
voids in the presence of primordial NG. More precisely, we are after the differential comoving number density 
dn com / dR com of voids of comoving size R com which satisfies 
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where R is the Lagrangian radius of the void related to the comoving radius by R com = 1.7 R (see below), 
S = <y 2 {R) is the variance of the density contrast smoothed on scale R, and the multiplicity function f(S) is 
what we will calculate analytically based on the excursion set formalism. Our results will be in a form that 
should be testable in A-body simulations which have complete control on the dark matter distribution. We 
will leave the second (and important) part of the exercise, namely that of relating our results to observationally 
relevant quantities, to future work. 

We should mention that this approach has been adopted by other authors before us [47, 48] (see also 
Ref. [49]). The analysis of Kamionkowski et al. [47], while pioneering, was based on several simplifying 
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assumptions which we believe give an incomplete picture of the problem. In particular they used a Press- 
Schechter-like approach for the statistics which can be shown to miss certain scale dependent terms in the 
multiplicity / that arise from somewhat complex multi-scale correlations [50, 51]. Additionally their treatment 
of the non-Gaussianities involved an Edgeworth-like expansion (based on Ref. [28]) for the one point probabil- 
ity density of the smoothed density contrast, linearized in the non-Gaussianity amplitude. As pointed out in 
Ref. [51] for the case of dark matter halos, such an expansion tends to underestimate the effect of the NG on 
large scales, which is the regime we wish to probe using voids. And finally, Kamionkowski et al. also ignored 
a complication which arises when studying voids, which Sheth & van de Weygaert termed the "void-in-cloud" 
issue (which we discuss in detail below) . Lam et al. [48] on the other hand used a more rigorous approach which 
accounted for multi-scale correlations and the void-in-cloud problem, but their treatment was also based on a 
linearized Edgeworth-like expansion and is hence subject to the same caveat mentioned above. Additionally 
their analysis made a technical assumption regarding the multi-scale correlations which is not strictly valid 
(see below). Our treatment below will be based on techniques introduced by Maggiore & Riotto [52, 53] and 
developed further by D'Amico et al. [51]. This will firstly allow us to derive a void multiplicity function which 
is valid on larger scales than the expressions of other authors. Secondly our approach will allow us to carefully 
account for the void-in-cloud issue (without making assumptions regarding the multi-scale correlations), and 
we will show that its final effects can actually be described using a simple yet accurate approximation. 

The plan of this paper is as follows. In Section 2 we review the excursion set analysis of Sheth & van de 
Weygaert based on the spherical ansatz for Gaussian initial conditions, and describe the void-in-cloud issue 
mentioned above. To carefully account for this problem, we turn to a path integral description. We begin in 
Section 3 by rederiving the Gaussian result for the void multiplicity using path integrals, which allows us to 
introduce some formal machinery as well as discuss some of the subtleties of the calculation in a controlled 
setting. Following this, in Section 4 we generalize the calculation to the non-Gaussian case and derive our main 
result for the non-Gaussian void multiplicity, showing how to account for the void-in-cloud issue. Our final 
result is given in Eqn. (68) and is plotted in Figs. 2 and 3. We conclude in Section 5 with a brief discussion of 
the result and prospects for future work. Several technical details have been relegated to the Appendices. 

2 Spherical void statistics: the SvdW result 

The spherical ansatz for an expanding underdense region allows for a completely analytical treatment, exactly 
like what happens in the usual spherical collapse model [41]. Detailed numerical treatments [54, 33] also indicate 
that, although simplistic, this ansatz goes a long way in describing the evolution of individual underdense 
regions in realistic settings. For mathematical details of the model we refer the reader to Ref. [37], where 
Sheth and van de Weygaert (henceforth SvdW) also extensively discussed the physics of voids in the context 
of hierarchical structure formation. Since our primary concern in this paper is the statistical model built using 
the spherical ansatz, we will restrict ourselves to briefly describing some of the physical aspects relevant to our 
treatment. 

An important feature of the spherical expansion scenario is the phenomenon of shell crossing [36, 37]. 
Physically shell crossing occurs because inner mass shells decelerate slower than the outer ones. With a 
sufficiently steep initial profile, this difference in decelerations is large enough that the inner shells catch up 
with the outer shells in a finite time, piling up mass in a very sharp ridge-like feature which then evolves 
sclf-similarly [43]. It can then be argued that the time of the first shell crossing is a sensible definition of 
the time of void formation. Such a definition brings us square into the realm of barriers and random walks, 
with the linearly extrapolated overdensity <5 C from spherical collapse being replaced by a linearly extrapolated 
underdensity (— <5 V ) corresponding to the time of first shell crossing. For an Einstein-de Sitter cosmology and 
assuming a tophat initial profile for the underdensity, this shell crossing threshold evaluates to 5 V = 2.72 [36]. 
Note that the excursion set ansatz relies on smoothing the (linearly extrapolated) initial conditions with a 
Lagrangian scale R associated with a conserved mass M oc R 3 . The observationally relevant scale however is 
the comoving scale R com which corresponds to the physical size today of a sphere reaching first shell crossing. 
The spherical model predicts the relation [36] 

-Rcom = 1-7R . (2) 
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It might seem that the statistical problem involving voids is now identical to the one with collapsed halos, and 
that a simple replacement S c — > — 5 V in all halo statistics results would suffice to give the corresponding void 
statistics. The situation is not this simple however, as SvdW discuss extensively in Ref. [37]. Recall that the 
collapse problem had to deal with the subtlety of the "cloud-in-cloud" issue. This was the fact that, when 
counting the fraction of collapsed objects of mass M, one must only consider those trajectories (of the random 
walk of the smoothed density contrast) which first crossed the threshold <5 C at a smoothing scale M, as the 
smoothing radius is decreased. Physically this corresponds to excluding overdense regions which are embedded 
in bigger overdense regions, since in such a case only the bigger region would survive as an independent collapsed 
object. Due to the nature of the problem, of course a similar requirement also holds in the case of voids, which 
SvdW call the "void-in-void" issue. This can be dealt with exactly as in the collapse case, and requires us to 
consider only first-crossing scales. 

Additionally, SvdW also identified a second problem which is unique to the case of voids. Physically, this 
is the scenario where a region of size i?i satisfies the threshold underdensity requirement and is a potential 
void candidate, but happens to be embedded in an overdense region of size R2 > Ri that satisfies the collapse 
criterion. This bigger region will then form a collapsed object, crushing the undcrdense void candidate out 
of existence. Explicit examples of such cases in TV-body simulations were shown in Ref. [37]. The problem is 
therefore to exclude such situations from the statistics and was labeled the "void-in-cloud" issue by SvdW, 
who also showed how the problem could be tackled. The basic idea is that for the void formation problem 
there are now two barriers which are relevant - the negative void shell crossing threshold (— 5 V ) which we will 
call the "void barrier", and also the positive halo formation threshold S c which we will call the "halo barrier" 1 . 
Statistically, excluding the void-in-cloud cases amounts to counting only those trajectories in which the void 
barrier is first crossed before the halo barrier is ever crossed, as the smoothing radius is decreased from large 
values. SvdW showed that accounting for the void-in-cloud issue qualitatively changes the behaviour of the 
void multiplicity as compared to the halo multiplicity, by introducing a cutoff at small smoothing radii. This 
is again intuitively clear since small underdense regions would be more likely to find themselves embedded 
in larger scale overdensities, and should therefore not survive as voids. Our interest when extending the 
problem to the case of non-Gaussian initial conditions will primarily be in the large radius end, which is where 
non-Gaussianities are expected to play a significant role. 

Before proceeding, a word on notation. In addition to the void barrier (— 5 V ) and the halo barrier 5 C , we 
will frequently also require: the total barrier height 

S T = S C + S V . (3) 

Additionally all our later expressions will involve the "dimensionless" parameters 6/S 1 ^ 2 , where S = S(R) = 
(S R ) is the variance of the density contrast smoothed on scale R. We will use the notation 

u c = 6 c /S^ 2 ; i, v = (5 v /S 1/2 ! ^t = <5t/S 1/2 . (4) 

In the following we will frequently refer to the variance S as "time" , in keeping with standard terminology in 
excursion set theory. Statistically we are therefore interested in the probability distributions of first crossing 
(f.c.) times, or f.c. rates (conditional or otherwise), which will be denoted by F with appropriate subscripts. 
The multiplicity / which appears in the mass function (1) is related to T by / = 1ST . SvdW showed that the 
conditional f.c. rate of the void barrier, accounting for the void-in-cloud issue and assuming Gaussian initial 
conditions, is given by 

00 

^Svdw(S) = E 3" sin CMv/*r) e-^ s '^ , (5) 

j=l °T 

whose Laplace transform is 

r°° Q sinh (V2sS c ) 

C Svdw (s) = / dSe- sS F Svdw (S) = . ; j 7=7^ • (6) 
Jo sinn (v2sdTj 

1 The actual value of S c must be chosen carefully, since simply setting it equal to the usual spherical collapse value 
of 1.686 physically amounts to only excluding those underdense regions that would be completely crushed by collapsing 
overdensities. As SvdW discuss, this all-or-nothing approach fails to account for intermediate underdensities which are 
in the process of being crushed, thereby overestimating e.g. the typical void size in the Gaussian case (see also Ref. [46]). 
In this work we will ignore this complication and assume 5 C = 1.686. 
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SvdW derived this result based on probabilistic arguments and a clever use of Laplace transforms. This result 
is also known in the condensed matter literature on first crossing problems, and can be found derived e.g. in 
Ref. [55], based on a solution of the Fokker-Planck equation 

d s n = ±%n, (7) 

where IT(<5, S) is the probability density for a diffusing particle in the presence of absorbing barriers at 6 = S c 
and 5 = —S v . 

In the following we will be mainly interested in the large mass or small S limit of Tsvdw- As it stands, the 
expression (5) is not particularly useful in this limit, since for small S an increasing number of terms become 
important 2 . It is very useful therefore to recast this expression in a form that is better behaved as S — >■ 0. 
This is not hard to do, and we show in Appendix B.l that the conditional f.c. rate .Fsvdw can be written as 

oo 

•Fsvdw(S) = — _ J2 K - 2j*r)e-<*- W 2 . (8) 

^ ' j=— oo 

(One can check that the series above is also identical to the one in Eqn. 46 of Lam et al. [48], with our j = 1 
term corresponding to their n = 1, our j = —1 to their n — 2, and so on.) As S — > 0, the j = term rapidly 
becomes the most important and the resulting multiplicity reduces to 

S^O : /svdwK^T) = 25J- Svdw (5)^^ v e-^ 2 /2 = /ps(l/v)j (9) 

which is just the 1-barrier Press-Schechter result for the void barrier. This is not surprising, since for small times 
it becomes increasingly unlikely for a trajectory to cross the halo barrier and return to cross the void barrier, 
reducing the result to a single barrier one. Figure 7 of Ref. [37] illustrates this effect, with the void-in-cloud 
effects becoming significant only at u v < 1.5 or so. 

In our non-Gaussian extension we would like to address the question of whether this intuition continues to 
remain true, or whether the void-in-cloud issue is now relevant at very large radii. If we simply assume that 
the Gaussian reasoning still holds, then the large radius end should be describable as a single barrier problem. 
This was essentially the reasoning of Kamionkowski et al. [47], who applied the Press-Schechter approach to 
a single void barrier. If we follow this reasoning, then at the least the non-Gaussian void multiplicity / vo ids 
should incorporate the single barrier effects discussed in Refs. [50, 51], and we should expect that / vo ids is given 
by replacing S c — > — S v in the expression derived by D'Amico et al. [51] (their Eqn. 68 with k — and a = 1 
to remain within the sharp-k filter and fixed barrier approximations we are using here). 

In this paper we will address this issue rigorously using path integral techniques. We will see that, while 
the void-in-cloud issue introduces some technical complications, the end result is in fact that at sufficiently 
large radii one can treat the problem using a single barrier. For smaller radii, we will see that the effects 
of the void-in-cloud issue are not negligible but can be described using a simple yet accurate approximation. 
Our results are broadly in agreement with those of Lam et al. [48] , although our final expression for the void 
multiplicity is different from theirs (see below for a more detailed comparison with their work). We begin in 
the next section by introducing some path integral machinery. 



3 Deriving the SvdW result from path integrals 

Our goal in this section is to reproduce using path integrals the SvdW result which assumes Gaussian initial 
conditions. This will allow us to rigorously extend the result to the non-Gaussian case. Unfortunately the 
calculations we will end up doing are technically rather involved. It is instructive therefore to first go through 
an intuitive derivation of J-"s v dw- 

Consider the probability distribution n a n owoc j (<5, S) for the location of the diffusing particle, which solves 
the Fokker-Planck equation (7) in the presence of two absorbing barriers and is nonzero in the "allowed" region 

2 The expression (5) is more useful in the large time (S — > oo) limit, which is typically encountered in condensed 
matter systems [55]. 
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— 5 V < 5 < 5 C between the barriers. One can interpret the rate T at time S as the amount of probability 
leaking across the chosen barrier per unit time. This is similar to constraining the probability to stay within 
the barriers until time S and then letting it evolve freely, as if the barriers were removed at time S. The rate 
at time S across let's say the void barrier can then be computed as a derivative of the probability that has 
leaked across this barrier at S, i.e. 

-<5v 

•Fir 



,(s) = d s f " dsu(s,s). (io) 

J — oo 



Here II represents the unconstrained probability density, which also satisfies the same Fokkcr-Planck equation 
(7), with the "initial" condition that at time S we have II = n a n wod- Using the Fokker-Planck equation one 
can simplify the expression above to find 



-^intuitive (S) — — dsHaUowed (S, S) 



(11) 

<5=-5 v 



The solution to the Fokker-Planck equation which vanishes at the two barriers S = S c and 8 = — S v and starts 
from a Dirac delta at the origin at the initial time is simply an infinite sum of Gaussians with shifted mean 
values, 




n a iiowcd(<5, S) = y~] < exp 



(8 + 2jS 



T J 



\2 



25 



exp 



(8 - 2jS T + 2S V ) 2 
2S 



(12) 



This is easy to verify by inspection, since each Gaussian is separately a solution of the Fokker-Plank equation 
(7), and by construction the infinite sum satisfies the boundary conditions Yl(—S v ,S) = H(S C ,S) = (as 
can be immediately verified for (— 6 V ), while for S c it is enough to shift j to j + 1 in the second Gaussian). 
Straightforward algebra then shows that the conditional f.c. rate computed according to Eqn. (11) is precisely 
the series given in Eqn. (8), 

^intuitive (S) = ^SvdW (£) . (13) 

These ideas can be made more rigorous using the language of path integrals. The path integral approach to 
excursion sets for computing the halo multiplicity was developed by Maggiore & Riotto (MR) in a series of 
papers [52, 53, 50], and was improved upon by D'Amico et al. [51] for the non-Gaussian case. In the following 
we will mainly refer to the techniques developed in Ref. [52] , restricting ourselves to the simplest case of a fixed 
barrier and a sharp-k filter (so that one is dealing with a Markovian stochastic process). The basic quantity 
one deals with is the probability distribution function W(8o; {(*>&}„; S) for a discrete random walk with n steps 
{^fc}|fe=i! where Sk denotes the matter density contrast smoothed on a scale Rk corresponding to a variance 
Sk , with steps of equal spacing AS in the variance starting at So with corresponding density contrast So (both 
of which we will assume to be zero) and with the last step denoted by S n = S. We have 

W(S ; {S k } n ; S) = ( 5 D (& - ft) . . . fo(4 -<*„)>, (14) 

which for the Gaussian, sharp-k filter case reduces to 

= J] * AS (4+i - S k ) ; * As(x ) = (2 7 rA5)- 1 / 2 e- a;2 /( 2AS ) , (15) 

fc=0 

with the superscript standing for "Gaussian, Markovian". Throughout our calculations we will deal with 
the path integral of this probability density over the "allowed" region - in the single barrier case that MR 
considered this would be the region Sk < S c , 1 < k < n — 1, while in our two barrier problem it will be the 
region — S v < Sk < S c . We will therefore be interested in objects of the type 

U AS (S ,5 n ;S)= [ dS^.-dS^WiSo-ASkUS). (16) 

J allowed 

This object is the probability density for the diffusing "particle" to remain inside its allowed region until time 
S. Our goal will be to calculate the rate at which this probability leaks out of the allowed region (i.e. the 
rate at which the "particles" escape) across one of the boundaries. Making these ideas rigorous requires us to 
introduce some technical aspects of the path integrals, and it will be easiest to do this in the more familiar 
single barrier case. Let us therefore begin with a brief recap of the (Gaussian) MR calculation for halos in the 
case of a fixed barrier and a sharp-k filter. 
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3.1 Recap of the single barrier problem 

As mentioned above, it is useful to define the constrained probability density n^ m bar as {5q, 5 n ; S) for the density 
contrast at step n, given that at all previous steps the walk remained below the barrier, so that 

n ? m bar,As(^oA;S)= dS, . . . dS n ^W^(6 ; {6 k } n ; S) . (17) 

J — OO 

MR showed [52] that the continuum limit of this quantity recovers the result of Bond et al. [40], 

n? m bQr , A ^ (^ -5; S) = u Bond (s , S;S)^^L= ( e -<«-*>) a /w _ e -(^-ss o) > (18) 

being the solution of the Fokker-Planck equation (7) with initial condition n(<5o, 8; 0) = Sd(S — So) and boundary 
conditions H(So, S c ', S) = = U(Sq,S — > — oo; S). We will go through the derivation of a similar result for the 
two barrier case below. 

We are looking for the distribution of the f.c. time S c at which the random walk first crosses the barrier 
S c . This can be constructed as follows. Start with the cumulative probability P(S C > S) that 5 C > S, which 
is the same as the probability that the barrier has not been crossed until time S, i.e. 

P(S C >S) = lim d,5„n? m bar (Jo, S n: S) . (19) 

In the continuum limit this quantity is straightforward to compute using Eqn. (18). The f.c. rate is then just 
the negative derivative of P(S C > S) w.r.t S, leading to the multiplicity 

/psK) = 2SJp S (S) = -2Sd s P(S c >S) = \j^v c e"^ 2 , (20) 

V 7T 

which is the celebrated Press-Schechter result (accounting for the infamous factor of 2), where v c was defined 
in Eqn. (4). 

On passing to the two barrier problem, in order to calculate the required constrained f.c. rate, the cumu- 
lative probability we need is P(S V > S,S C > S v ), which is the probability that the void barrier is first crossed 
before the halo barrier is ever crossed. In this object, the conditioning of the stochastic variable S v is thus 
on another stochastic variable S c which is the first crossing time for the halo barrier, and this complicates 
the issue. In the Gaussian case of course one can solve the problem the SvdW way, without resorting to path 
integrals. We are ultimately interested in the non-Gaussian case though, and we therefore explore a path 
integral derivation of the Gaussian SvdW result, which we can generalize later to the non-Gaussian case. 



3.2 ^svdw from a path integral analysis 

The basic trick we employ is to exploit the discretized nature of the path integral, and to pass to the con- 
tinuum limit carefully. To keep the discussion as clear as possible, we relegate several technical derivations 
to the Appendices. As in the single barrier case, it is useful to construct the constrained probability density 
n^(<5o, S n ; S n ) for the density contrast at time step n, given that the walk has not crossed either barrier at 
any previous time step, 

n s ™ s (5 ,S n ; S n ) = J 5 ^ dS 1 . ..4-^ gm (^o; {**}„; S n ) . (21) 

We will need to compute this quantity and its integrals under various limits and assumptions, which we will 
come to presently. To begin with, note that with a discretized variance parameter S one can explicitly construct 
a probability (rather than a density) for the constrained first crossing of the void barrier to occur at a specific 
time step n. Namely, the integral J_ Sv dd n Tl^(do, S n ; S n ) is the probability that the walk did not cross either 
barrier for the first n — 1 steps, and crossed the void barrier at step n. This is the same as the probability that 
the first crossing of the void barrier is at step n, and that the halo barrier has not yet been crossed. Using 
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this, the cumulative probability P(S V > S,S C > S v ) is obtained by summing over all possible choices of the 
final step n that have S n > S, and then passing to the continuum limit, 

/— <5 V /•oo i p — <5 V 

dS n U A m s (S ,8 n ;S n )^ dS Urn — / dS n U A m s (5 ,d n ; S) . (22) 
.. , -oo JS AS^O A£ J_ oa 

The required f.c. rate is simply the negative time derivative of this object, which we can read off as 

Jws(S)= lim -L f ^ dS n U A m s (5o,5 n ;S), (23) 

which involves the integral of on the "wrong side" of the void barrier, 6 n < —5 V (hence the subscript WS). 
One can also see why this is the correct object to compute, since it corresponds to the amount of probability 
leaking out of the void barrier in a time interval AS, divided by AS. While the integral itself will vanish 
in the continuum limit, we will see that it does so like ~ AS, leaving a finite limit in Eqn. (23) such that 
Fws(S) = J^SvdwiS)- 

The above arguments may seem like a convoluted way of arriving at the result. When we move to non- 
Gaussian initial conditions for the double barrier however, we are not left with much choice in the matter. 
Nevertheless, it would be reassuring to know that the "wrong side counting" described above actually works 
in some other situation which is under better control. The single barrier calculation provides us with such a 
check. In this case, we can apply exactly the same arguments as above to obtain the f.c. rate, and we have an 
alternative derivation of the rate due to MR. As a check therefore, we should find that these two methods lead to 
the same answer. In Appendix B.2, we show that this is indeed the case for completely general initial conditions 
(i.e. without assuming Gaussianity) . In Appendix B.3 we also compare our approach with that of Lam et al. 
[48]. Their analysis effectively makes the assumption that the probability distribution W(So; {Sk} n ', S) defined 
in Eqn. (14) is factorisable, which is certainly true in the Gaussian case for the sharp-k filter (see Eqn. (15)) 
but is not valid e.g. in the presence of non-Gaussianities. We show in Appendix B.3 that our expression for the 
rate (both the Gaussian one of Eqn. (23) as well as its non-Gaussian generalisation which we discuss in Section 
4 below) is identical to what they aim to calculate in Eqn. 41 of Ref. [48], but without making this assumption 
of factorisability. With these reassurances, we proceed to the main two barrier calculation of this section. Since 
the calculation is rather technical in nature, the reader who is willing to take the result Jws = •T'Svdw on faith 
may skip directly to the non-Gaussian generalization in Section 4. 

To calculate the integral in Eqn. (23), we start by exploiting the factorisability of W sm to write 

U A m s (6 , 6;S + AS) = dx^ AS (x)U A m s (S , 6-x;S) 

Js-s c 

\ L ) an 



= E 



n=0 



E 

n=0 



d^U A m s (S ,S;S) f ' dx9 AS (x)x n 
Js-s c 

^—L d ^ s (6 ,5-,S)1^ S (5), (24) 



where the first equality follows from the definitions of H A ™s, W sm and some relabeling of dummy variables, the 
second follows from a Taylor expansion and exchanging the orders of integration and summation, and the last 
line defines the functions T^ s {8), which reduce to 

X^{S) = ^-f dye-* 2 y*. (25) 



(5-5 c )/(2AS)i/2 



Eqn. (24) is completely general, and holds for arbitrary 5. The first thing we can check is that this expression 
implies that in the continuum limit, n^ =0 (#o, S; S) vanishes on both barriers 5 = — 5 V and S = 5 C , 



nS=o(*o, Sv, S) = = n^ =0 (<5 , S c ; S) . (26) 

This follows by Taylor expanding the l.h.s. of (24) for small AS and comparing the lowest order terms on both 
sides, using 1 A S ^°(S C ) = 1/2 = l^ s ^°(— 5 V ). Here we set 8 to S c or (— S v ) before taking the continuum limit. 
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Figure 1: Schematic view of the two barriers (vertical dashed lines) with corresponding boundary layers 
(vertical dotted lines). The boundary layer functions v(rj) and 7/(77) are also shown schematically. 



This is a subtle point since these operations do not commute, which is clear from the structure of 2^ s (6). This 
brings us to the issue of "boundary layer" effects. 

In evaluating the integral in Eqn. (23), we wish to take the continuum limit after integrating, in contrast 
to what we did above for I$ s . The integral therefore depends not only on the behaviour of far from the 
void barrier, but also on its detailed behaviour in the "boundary layer" where \S + S v \ ~ %/ AS. In our wrong 
side counting approach, there are actually two boundary layers we must worry about. If we define the quantity 
V by 

(27) 

V2AS 

then the two boundary layers correspond to regions where \rj\ ~ 0(1), with the one discussed above corre- 
sponding to rj < 0. We will also deal with the second layer at rj > 0, which is in fact similar to the one discussed 
by MR in their single barrier calculation. In principle there are two more boundary layers, on either side of 
the halo barrier as well, but we will only need to deal with the ones near the void barrier. Fig. 1 illustrates the 
situation. To account for boundary layer effects in the integral in Eqn. (23), we use the approach discussed 
by MR in Ref. [52]. This involves two steps: first, we calculate the leading behaviour of n^(<5o, 8; S) for 
small AS when we hold 6 < —5 V to be fixed - i.e. in the limit of large negative r\. Denote this function by 
C^g(So,S;S). Next, in order to calculate ±1^, for fixed but small AS* and arbitrary S < ~S V (so that 77 is 
negative and arbitrary), we introduce a boundary layer function v{rf) which we choose to normalize so that 
v(rj — > —00) — > 1, and then write 

n^(<5 , 5; S) = »(»,)cg {So, 5;S) + ..., (28) 

where the ellipsis indicates terms of higher order in AS* than the leading order in C^. . The function v(rf) will 
in general behave nontrivially in the boundary layer where |?7| < C(l), capturing effects missed by the function 

^AS on ^ s own - 

Let us now calculate C^g which can be done by starting with Eqn. (24). In Appendix C.l, we show that 
in this limit of large negative rj the functions (5) reduce to 

1/2 



2^ < -m = (f ) ' (« + «v)" -1 (1+0 (^y ) 



(29) 



where the ellipsis indicates terms which are exponentially suppressed. Replacing this in Eqn. (24) we find that 
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to leading order in AS, the summation over n can be carried out exactly using the identity 



£ ^J- {6 + S v ) n ^n| m s (<5 , S; S) = UTsiSo, S v ; S) , (30) 

n=0 

which is not zero for finite AS*. To leading order in AS, the l.h.s. of Eqn. (24) is simply the function 
C^g(do,S; S) we are looking for, and we can set 

/ASA 1/2 P -(<5+M 2 /2AS 
c£>(5o,5-,S)={—) ' ls + §vl UZ(5 ,-5 V ;S), (31) 

where it is understood that we are only interested in the leading order in AS. Since the object n^(<5 , — S v ; S) 
is independent of 5 we will calculate it later. The function II^(<$o, 6; S) for arbitrary S < —5 V can be written 

as 

nT s (So,S < ~S V ; S) = 1 "^ e-^IlS^, -<5 V ; S) , (32) 

where we have expressed C^g in terms of rj. In the other limit of fixed AS but 5 — > — <5 V from below, due to 
continuity of 11^ (5q, <5; 5*) we have 

v(rj) -> 2^(-ri) as 77 -> 0" . (33) 

Let us focus now on the integral in Eqn. (23). It is not hard to see that this can be re-expressed as an integral 
over rj from —00 to 0, and reduces to 



^n r ^o^;S) = (-j 



w(?y)e ' 



nS(io,-i5,;S). (34) 



The integral in square brackets depends strongly on details of the function v(n), and is in fact dominated by 
contributions from the boundary layer where \n\ < 1 due to the exponential cutoff. We do not know of any 
way to calculate this object from first principles. However, this integral is a finite constant since the integrand 
is well behaved as n — > (see Eqn. (33)) and is exponentially cutoff as n — > —00. Additionally, this constant 
only depends on the properties of the single void barrier. We will see below that it can then be fixed by using 
our knowledge of the single barrier calculation. 

The calculation of Jws wn l be complete once we evaluate the object lT^,(<5o, — S v ; S) (which we need only 
to leading order in AS). We start by returning to Eqn. (24), where we are now interested in the limit S — > — S v 
from above at fixed AS - i.e as n — > + . As we shall see next, this will allow us to exploit the properties of 
the continuum probability density n^_ which is non-zero in — S v < S < S c and is moreover calculable. As 
one might expect, once again we need to account for a boundary layer, this time the one where rj is positive 
and of order unity. This calculation closely mimics that of a similar object in the single barrier case [52] . We 
essentially repeat the previous analysis of Il^(do, 6; S), now on the other side of the void barrier. Once again 
we define a function C^g (So, S; S) which reproduces the leading behaviour in AS of LT^((5o,(5;S) when we 
hold S e (— S V ,S C ) fixed and let AS — > 0. This time we can easily see that C^g (So, S; S) must simply be the 
continuum limit probability density 

C^(S ,S;S)^n^o(So,S;S). (35) 

Introducing a second boundary layer function u(n), we can express II^(<5o,£;S) for arbitrary S € [— S V ,S C ) as 

IlT^o, 5; S) = u{ri)C${6 , 5; S) . (36) 

With C^g(So, 5; S) defined to be smooth in the boundary layer where n = 0(1), we can obtain its leading 
behaviour in AS by a simple Taylor expansion, which leads to 

(So, S; S) = n| m s=0 (<5 , -S v + nV2AS; S) 

= rjV2ASd s U^ s=0 (6o,S^-8 v ;S) + ... , (37) 
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where we used the fact that the continuum limit function n^ =0 vanishes when evaluated at the barrier. 
Finally, we can write the leading order form of 11^ (Jo, —5 V ;S) by taking the limit n — > + in Eqn. (36) to get 



H^ s (So,-S v ;S)^ 1 V2ASd s nZ =0 (So,S = -S v :S) ; 1 = limqufa). (38) 

r;— >0+ 

Here 7 is another constant which depends on details of the boundary layer near the void barrier, and is in 
fact very similar to the constant which MR compute in Rcf. [52] (see their Eqn 78). We will however simply 
combine it with the other unknown constant which appears in Eqn. (34) and fix it below by appealing to the 
single barrier problem. Plugging the above relation into Eqn. (34), we see that the integral in Eqn. (23) is 
proportional to AS at the leading order. Putting everything together, the two barrier conditional f.c. rate is 
then given by 



T ws =A^d s U^(S ,S = -S v ;S) ; A = / -^-v(n)e-^ . 

1 V 71 " J -co {-V) 



(39) 



We now need the derivative of the continuum limit solution between the barriers. If we hold 6 € (— d v , S c ) 
fixed and take the limit AS — > then it is easy to see, by expanding both sides of Eqn. (24) in powers of AS 
and equating the lowest order terms, that the continuum limit function satisfies the Fokker-Planck equation 
(7) with boundary conditions (26) and initial condition a Dirac delta centered at So (which we set to zero for 
convenience). In other words, our rigorous derivation has reproduced Eqn. (11), upto a constant factor. Since 
we saw earlier that the expression (11) is the same as Jsvdw m Eqns. (8) and (5), we have 

T WS (S) = AT SvdW (S) . (40) 

We now show that the constant A must be unity. It is clear from the derivation above that A depends only 
on the properties of the void barrier. In particular, sending the halo barrier to infinity, S c — > 00, would not 
change the value of A. This, however, is precisely the limit in which one is solving the single barrier problem. 
We already know from Appendix B.2 that the single barrier wrong side counting argument for the f.c. rate is 
equivalent to the MR derivation, i.e. J-"ws,i-bar = J-MR.i-bar- Additionally, a calculation identical to the one 
above (but in the limit of S c — > 00) shows that .Fws,i-bar = -4.?MR,i-bar, and hence we find A = 1. 

4 Non- Gaussian voids 

With the path integral machinery in place, it is formally staightforward (although still somewhat involved in 
practice) to extend the calculation to the case of non-Gaussian initial conditions. We simply replace W &m and 
Has from the Gaussian calculation with their appropriately generalized versions. Namely, we have [50] 

W(S ; {S k } n ; S) = ( -Si)... 6 u (S n - S n ) ) 
= exp 



1 . . . 

C 5 3 5 kk)cdjd k di 

' j,k,i=i 

1 ™ .... 
+ 4! X! ( S :i s k$i$Tn )cd j d k d l d m + ... 



j,k,l,m—l 



W^(S ;{S k } n ,S), (41) 



where ( SjSkSi ) c , (SjSkSiS m ) c , etc. are the connected moments correlating different length scales (which all 
vanish in the Gaussian case, giving back W gra ). Using this, we also have 



/ dS-i . . .dS n -iW(S ;{S k } n ;S n ) . (42) 

One can now apply exactly the same wrong side counting arguments as in the Gaussian case (see the application 
to the single barrier case in Appendix B.2), and find that the required non-Gaussian f.c. rate is given by 

Jn G (S)= lim -L / V dS n U AS (6 ,5 n ;S). (43) 
as^o AS J_ 00 
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We will focus on weak primordial non-Gaussianities (NG) which are expected to be generated in inflationary 
models. One expects that in this case a scale dependent perturbative treatment of the NG along the lines 
discussed by D'Amico et al. [51] for the halo multiplicity should work well for voids as well, with possible 
complications due to the void-in-cloud issue. We proceed as in Ref. [51], beginning by defining the "equal 
time" correlation functions 3 (recall S = <rfj), 

£ „_ 2 = i% . n > 3 , (44) 

which are expected to follow the perturbative hierarchy s n ~ e™ for n > 1, where e <C 1, in generic inflationary 
models (see Appendix A for a brief introduction to primordial NG). Next we Taylor expand the "unequal time" 
correlators appearing in Eqn. (41) around the final time S [52] as e.g. 

(W>c= E Ja\r\ St 9,r) (S)(S-S j nS-S k y(S-S l y, (45) 

p,q,r=0 



where we introduced the (scale dependent) coefficients 



g£' q ' r) (S) 



d p d g d r 

d^dsidsf {6{Sj)s{Sk)s{Sl) ^ 



(46) 

- Sj—Sk—Si—S 



and similarly for higher point correlations. Such an expansion works well when considering large scales for 
which S is small. In this work we will restrict ourselves to the effects of at most C/3' , ignoring the higher 
order correlations which will be suppressed by powers of S and can be tracked as part of the theoretical error 
we make in the calculation [51]. We also introduce a convenient parametrization of Q^' ' ^ as in Ref. [51] by 
defining the function c\ (S) via 



gi 1 '°'°\S) = le 1 (S)c 1 (S)S 1 / 2 , (47) 
which satisfies the identity 

, „ ., 2 dhi£i , . 

Cl(5) = 1 + 3.hT<r (48) 

We will further follow the analysis of Ref. [51] and simplify the expression for the f.c. rate (43) by retaining 
equal time correlators in the exponential while linearizing the unequal time contributions, which is a strategy 
that works well at least for the single barrier problem. We then find, to the order we are interested in, 



Jng(S') = lim -!— f dS n [ dSn-i ■ ■ - ddiexp 

AS-s-0 A5 J_ OD J_ s 



1 , a, . ^ 1 



-{51)^9^ + -{5 4 n ) c £ 5A^ m + ... 



1 + lgi lfi '°\S) J2(S Sj)dj dkdi + . • . I . (49) 

In the single barrier case, the analysis at this stage was simplified by the existence of a very useful identity 
which states that, for any function g(5i, . . . , S n ), 



I 



s c n d / " 5, 



dSi... dS n ^2 d o9 = gf I d5i... dS n g . (50) 



3=1 



This identity is quite easy to prove, in fact it takes only a little thought to see why it should be true. Remarkably, 
there is a similar identity which turns out to be very useful for the two barrier case also, which is not as obvious 
to write down as (50). As we show in Appendix D.l, for any function g(d\, . . . , 5 n ) it is also true that 



dS n dS n --i . . .dSi y^djg = - — / dS n dS n -i . . . d5\g . (51) 

-00 J-5 V „ , a °v St J -00 J-5 V 



3 = 1 



3 In terms of the reduced cumulants <Sa, 54, etc. we have £1 = <r<S 3 , £2 = cr 2 <S4, and so on. 
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Notice that the final derivative holds fixed the sum St = S v + S c rather than S c alone. For brevity, throughout 
the rest of the paper we will omit the explicit reference to this fact, and simply write dg v in place of 9/<9<5 v |<5 T . 

We can immediately see why this relation is useful. Consider the equal time exponentiated derivative 
operator in Eqn. (49). The identity (51) allows us to pull this entire operator outside the integral, exactly like 
in the single barrier case. The remaining terms involving the integrals and the continuum limit can now be 
treated individually. We can recognize the first of these as simply J-"s v dWj an d defining the second as 

J-( 3 > NL > ee A lhn o J ^ dS n d8 n ^ . . . d5i\g£ fi '°\S) £(S - S^dj ]T M W*» (52) 

3 k-J 

(with the notation (3,NL) standing for 3 point, next to leading), the result at this order is 

J NG = e (l/3!)e 1 S^ai + (l/4 ! ) e2 5^i+- (jr^ + J-(3,NL) + _ _ \ (53) 

The term J"( 3 . NL ) is tricky to evaluate, and in Appendix D.2 we show that it reduces to 



^ 3 - NL > = -^—e lClS ^dl f S -^L^ SvdW (S) 
4V2^ Jo ^S-S 



(54) 



where we wrote CJg 1 ' 0,0 ' in terms of c\. Note that .Fsvdw does depend on both <5 V and St, even though this is 
hidden by our compact notation. From Eqn. (8) we can write 

^Svdw(S) = -=- V B,e- B V 2 ; B j = u v -2jtA r , (55) 

V 27T O — 

] = -co 

and the integral over S in Eqn. (54) can be evaluated exactly, giving 

f ^^Svdw(S) = f] -^r 2 e- B > /2 , (56) 
Jo ^S-S j^oo S/ 

using which we get an expression for J"( 3 > NL ) ; 

J-< 3 < NL > = -_^i £lCl y (B]-\)er B V\ (57) 



Further, noting that in Eqn. (53) the combination S x l 2 ds w can be written as d/du v \ VTi s = d Vv , we can write 
J"ng as 

-F NG = 1 e(V30e 1 a3 v+ ( 1 /4,)e 2 a^ + ... g c - fl »/ a U 1^ {b] l} + _ \ (5g) 

v j=— oo ^ ' 

Consider the j = term, for which B — v v . The exponential derivative can be computed using the saddle 
point approximation as discussed in Ref. [51], and gives precisely the single barrier f.c. rate computed there 
for the fixed barrier excluding filter effects (sec their Eqn. 60 with v — > —v v ), 



2 e (l/3!)£i9j v +(l/4!)e 2 9* v +. 



("v - iyl ~ 1) + • • 

= y|^ e -^v 2 (i+ £l -v/3+(,?-, 2 /3), v V4) ^ _ l £iI/y(ci _ 4) + . . . 

= /NG,l-bar(fv) • (59) 

As discussed in detail in Ref. [51], the ellipsis in the second line denotes all terms that are parametrically 
smaller than the ones written down, which correspond to terms of order 0{ev~ x , e 2 ^, e 3 ^), where e -C 1 is the 
parameter controlling the NG via £ n ~ 6 n . We will assume that the largest scales we access satisfy ez/^ ~ 
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One can check that for /nl ~ 100 this corresponds to Lagrangian scales R <~ 25ft~ 1 Mpc, or very large voids 
(see Eqn. (2)). In this case one finds that the ignored terms listed above are all of the same order of magnitude. 
At all smaller scales, the term of the form gives the largest theoretical error. We will use this in our 

arguments below. 

For j 7^ 0, notice first that the Bj are all linear in v v , so that the effect of a derivative <9„ v |„ T is identical 



to that of a derivative w.r.t B y We therefore have 



i \ - e (l/3!)e 1 9|.+(l/4!) £2 a|.+.. 



j = -oo 



(60) 



At least formally, one might say that the solution is therefore just a series of single barrier results, leading to 
a multiplicity /ng = 25Jng given by 



/ngK^t) = J- V Bj e-^+^+^-^)B^ + ...) U_ 1 lB . (ci _ 4 ) + . 

V 7T z — ' \ 4 

j=-oo 

oo 

= /NG,l-bar(-Bj) , (61) 

J=— oo 

with = i^v — 2jVt and /NG,i-bar defined in Eqn. (59). The problem with this expression is that for any fixed 
v v , i>t and NG parameters £\, £2, etc., for large enough j the terms being ignored in the ellipsis will become 
comparable to the ones being retained. This is not such an important issue for the polynomial NG terms, which 

1 p2 

on their own would always be suppressed by the Gaussian factor j, This is also the reason why analyses 
such as those of Lam et al. [48] and Kamionkowski et al. [47] , which are based on the Edgeworth expansion 
and therefore have multiplicities of the form e~ v / 2 multiplying a polynomial in u, are not susceptible to the 
void-in-cloud issue for large enough voids. As D'Amico et al. [51] argued however, the Edgeworth series results 
break down when the combination eis 3 becomes of order unity. The D'Amico et al. analysis, which we have 
used here, instead effectively resums potentially troublesome terms and leads to the non-trivial series in the 
exponential in /NG,i-bar(^), which can give significantly different results for the halo multiplicity than the 
Edgeworth-like analysis (see e.g. Fig. 5 of Ref. [51]). 

In the two barrier case as well, these exponentiated terms are expected to be important for the j = piece 
at large enough u v . However for j 7^ they are problematic, being a series in the supposedly small parameter 
eBj. Clearly \eBj\ becomes larger than unity for large enough j (positive or negative) at any fixed v v , and the 
series expansions for all j values beyond this point break down. This is not surprising, considering that this 
expression depends on the saddle point approximation, which for each j is only valid provided \eBj\ < 1 [51]. 
A more intuitive way of understanding this breakdown is to note that the result for the f.c. rate is effectively 
a series of single barrier f.c. rates with successively larger barrier heights. Since each single barrier f.c. rate 
involves the behaviour of a non-Gaussian conditional p.d.f. (IIi_bar,AS=o) at the barrier, for successively 
larger j we are effectively sampling further and further extremes of the non-Gaussian tails of the distributions, 
which eventually can no longer be described perturbatively. In practice of course, we don't expect that the 
resummation of this series would dominate the Gaussian suppression to the extent of giving order unity features 
in the single barrier f.c. rate. 

In fact one can make a stronger statement based on the limit in which we artificially send S c — > 00, 
which must recover the single barrier result for the void barrier. In the expression (61), sending 5 C — > 00 
for fixed <5 V and S is the same as sending — > 00 at fixed v v . In this limit, if the single barrier result is 
to be recovered, then each term with j ^ must individually vanish. For a given j, as — > 00 we have 
Bj — > ±00 with the sign depending on the sign of j. A given term with j =/= will then vanish only if the series 
l + eiBj/3+ (e? - e 2 /3) B]/4 + . . . which appears in the exponential, resums as \Bj\ — > 00 into a form which is 
bounded both above and below by strictly positive numbers. In other words, the exponential suppression of the 
single barrier multiplicity must qualitatively remain intact. We will make the mild assumption that the lower 
bound on the resummed series is not arbitrarily close to zero but is closer to unity, which excludes pathological 
features such as e.g. sharp repeating spikes in the single barrier multiplicity with a slowly decreasing maximum 
height. 
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With this we can extend the argument to finite ut, by noting that in /NG,i-bar(-Bj) sending vt — > oo 
for fixed j is the same as sending \j\ — > oo for fixed ur- The previous paragraph immediately implies that 
the contribution of terms with increasing \j\ is progressively suppressed. Consider first a term in which the 
series expansion has broken down, so that \eBj\ = a n > 1 (while ev v is still significantly less than unity). The 
arguments above suggest that the quantity /NG,i-bar(-Bj) resums to the form 

hci-bUBj) ~ a 1 B j e-l a ' B ' = -a x aoe-^ a ^ , (62) 

where ol\ and a-i are positive numbers with a possible mild dependence on \Bj\- Let us compare this term with 
the biggest term crr we do not calculate in /NG,i-bar(-Bo), wrucn is 

err - ^-^(1+^/3+^/3)^/4) x Q ^-i) _ (63) 

The ratio of these terms is 

erro e 2 

and as long as this ratio is less than unity the error we make by ignoring /NG,i-bar(-Bj) is smaller than the 
one made by truncating the series in /NG,i-bar(-Bo)- We argued above that the quantity ot2 is strictly positive. 
Moreover we expect a\ to be of the same order as cx.2, since both are proxies for resummed series which are 
structurally similar. Since we cannot compute these objects, let us separately analyse the situations where a2 
and «i are of order unity or are much smaller. 

• If a.2, ai are of order unity, then since ao > 1 the first term in the exponential is much larger than the 
term containing (e^ v ) 2 , and clearly we have r<l. 

• Even if a 2 , Qi are not close to unity, in practice r remains small provided only that ct2 is not very close 
to zero, this assumption following from our comment above on the lower bound on the exponentiated 
series. For example if we assume 0:2, |ai| > (e^v) 2 then we find 

r < a ^e~^v(^-i) . ( 65 ) 

Now even if ao is not much larger than unity the ratio r will remain small for all interesting values of 
u v . Say ao ~ 1.5, then r is already less than unity at u v = 1, which is the lower limit for our formalism 
in any case. 

To summarize, j-values for which the series expansion in \eBj\ has formally broken down (even mildly) are not 
expected to give contributions larger than the terms which are already being ignored in the j — term. 

Finally, for terms in Eqn. (61) with non-zero j where the series expansion has not broken down, one might 
still expect that the contribution of these terms is always smaller than err . We will now see that this is not 
quite true, although it is possible to ignore the j ^ terms if we only probe the largest voids. We start by 
comparing the leading order piece of /NG,i-bar(-Bj) which is given by Bje^^ B ^ 1+ °^ B ^\ with err . The ratio 
of these terms is 



Bje 



-iBj(l+0( £ U,)) 



err 



3t e -*W-"*+-) , (66) 



with the ellipsis denoting parametrically smaller quantities. In order to be able to ignore this term, we need 
to place an upper bound on the above ratio and show that this bound is less than unity. Since this entire 
discussion is valid only if ,/nl (else we simply use /svdw), suppose now that we have a lower bound for |/nl|, 
say |/nl| > 1 in the local model, which translates to e > 2 • 1CP 4 or e -1 < 5000. Writing Bj = v v (1 — 2j8^/5 v ) 
with St/S v s=s 5/3 for S c ~ 1.7 and S v ~ 2.7, we can see that \Bj\ > 2v v for all j ^ 0. Together with the bound 
on e, this allows us to calculate a minimum value ^ vm in such that the ratio (66) is always less than unity. We 
do this by setting \Bj\ = 2v v since all larger values will give smaller ratios, and find that the ratio (66) is 
guaranteed to be less than unity if v v > ^ vm in — 2.5 corresponding to voids with Lagrangian radius larger than 
<~ 6.5/i _1 Mpc (or comoving radius larger than ~ ll/i _1 Mpc) in a IVM^lP-compatible ACDM cosmology [20]. 
If we assume a smaller lower bound on |/nl|, the value of ^ vm in wm increase slowly. 
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Figure 2: The ratio of non-Gaussian to Gaussian multiplicity for /nl = +100 (left panel) and /nl = —100 (right 
panel) in the local model, at redshift z = 0. The two barrier and single barrier ratios are indistinguishable, as discussed 
in the text. There is a clear departure at large radii from the ratio 7?kvj proposed by Kamionkowski et al. [47], which 
is the one calculated by Lo Verde et al. [28]. The error bars represent theoretical errors of order O^ev' 1 ) as discussed in 
the text. 



We do not need to be this conservative however. We can also consider values t-v < ^v m in by retaining 
an appropriate number of terms with j =^ 0, depending on the chosen value of /nl or e. In fact as we will 
see presently, in practice there is an even simpler way of accounting for the void-in-cloud effect. Firstly, our 
arguments above indicate that at a given v v > 1 we must keep the leading order behaviour (schematically 
~ Bje~2 B j( 1+eB i + --^) of those j ^ terms for which \eBj \ < 1 and the ratio in Eqn. (66) is larger than unity. 
(An analysis similar to the one above shows that the subleading terms can always be ignored in this case, 
so that the dominant error is still given by erro) This is straightforward to implement numerically, and in 
Fig. 2 we show the ratio /ng^vi v t)I /svdw^v, v t) (solid red) as a function of R com for /nl = ±100 at z — 
with S v — 2.72 and S c = 1.686. We also plot the ratio of the single barrier functions /NG,i-bar(^v)//ps( J/ v) 
(short-dashed black), see Eqn. (59). As we see, these curves are indistinguishable, and one can also check that 
even at the smallest radii we consider, the difference between the curves is less than 0.5%. In fact this result 
is also easy to see analytically, by recognizing that the j ^ terms give increasingly smaller contributions to 
the sum in /ng(^v, vt)- One then has 

/ng("v,Vt) _ /NG,l-bar(fv) + J2j^0 /NG,l-bar(-Bj) 



/svdwK, v T ) /psK) + Ej^o fps(Bj) 



/NG,l-bar 
/ps(^v) 



K) 



E 

j¥0 



/NG,l-bar(-Bj) fps(Bj 



/NG.l-bar(^v) /ps(^v) 



(67) 



where we linearized in the j ' ^ terms. Now, for large is v the summation in the second line of (67) will be 
suppressed simply because of the Gaussian factor e 2 i " v) in each term. For v v — > 1 this supression will not 
be very strong at least for small values of \ j\ ^ 0. However, for such terms the single barrier ratio /NG,i-bar//ps 
approaches a constant, leading to a cancellation of the terms in the summation above. In practice therefore, 
to extremely good accuracy the two barrier non-Gaussian multiplicity can be simply written as the product of 
the single barrier non-Gaussian ratio with the two barrier Gaussian multiplicity /svdw, 



oo 



/ngK,^) = y \ hG,l-bm(Bj) RJ ( '^ NG ' 1 ba '^ v ) ) x /svdw(^v,^T) , (68) 



j=—oo 



where /NG,i-bar(^) was defined in Eqn. (59), fpsi^) in Eqn. (20) and Bj = i/ v — 2jv^. Of course this discussion 
is subject to the caveat that there is always a theoretical error at least due to the terms we do not compute 
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Figure 3: The two barrier non-Gaussian multiplicity /ng as a function of u v , for /nl = ±100. The inset shows a 
zoomed in view of the range highlighted by the cyan (shaded) area, corresponding to the range R CO m = 2-4/i _1 Mpc. 
The curves in the inset correspond to /s v dw (gau), /^ ,L ~ ±10 °' ) (2±) and iZba, 00 ' (1±)- The single barrier result 
overestimates the void abundance by order ~ 10% for both signs of /nl- The correct, two barrier result on the other 
hand displays a behaviour opposite to that at large R CO m '■ positive /nl slightly enhances the abundance at small R CO m- 
Also shown are the theoretical errors on the two barrier result as discussed in the text. For clarity we only show these 
for ,/nl = +100. 

in the single barrier multiplicity. For comparison, in Fig. 2 we also show (dashed blue) the ratio proposed by 
Kamionkowski et al. [47], which (as expected) deviates from our prediction at large radii. The same will be 
true of the ratio calculated by Lam et al. [48] at large radii. Fig. 2 also illustrates the complementary nature of 
the void abundances as a probe of NG, since the abundance at large radii is reduced compared to the Gaussian 
case for positive values of /nl and vice versa for negative /nl , which is the opposite of what happens for halo 
abundances. 

In Fig. 3 we show the non-Gaussian multiplicity as a function of v v , together with a zoomed in view at 
small comoving radii (v v > 1). We see that ignoring the j ^ void-in-cloud terms entirely will overpredict the 
void abundance at small radii for either sign of /nl- On the other hand, accounting for void-in-cloud effects 
gives a result which is approximately the same as the Gaussian one, with a slight enhancement for positive 
/nl and a slight reduction for negative /nl- Notice that this is the opposite of what happens at large radii. 
Unfortunately the magnitude of the reversed effect at small radii appears to be too small to be observationally 
relevant. Finally, in Fig. 4 we show the differential comoving number density dn com /d log i? com defined in 
Eqn. (1) as a function of comoving radius i? C om- 

5 Discussion 

Primordial non-Gaussianity (NG) can be probed by the imprints it leaves on the late time large scale structure 
of the universe by modifying the distribution of matter, which is manifested for example in the abundance of 
collapsed objects. In this paper we have explored a second manifestation of this effect, which is the abundance 
of voids or underdense regions. While such calculations have been performed earlier in the literature [47, 48], 
they have been subject at least to the caveat that their treatment of the NG was based on a linearization in 
/nl, which is known to potentially misestimate the abundance of very massive objects [51]. Our calculation 
was based on path integral techniques introduced by Maggiore & Riotto [52, 50], and importantly also on the 
improvements to these techniques developed in Ref. [51]. The latter allow us to access a larger range of length 
scales than treatments based on linearizing in / N l, while the use of path integrals also allows us to carefully 
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Figure 4: Comoving number density denned in Eqn. (1), as a function of comoving radius i? CO m at redshift z — for 
/nl = ±100. 

account for the "void-in-cloud" issue first pointed out by Sheth & van de Weygaert [37], which is unique to 
the case of voids. We showed that in the final analysis, the complication introduced by the void-in-cloud issue 
can actually be accounted for in a fairly simple manner, with the void multipicity /ng^v^t) being given 
by the approximation in Eqn. (68). Our treatment of the multi-scale correlations arising from primordial NG 
is more rigorous than that of Lam et al. [48] (although these effects are small), and additionally our explicit 
expression for the multiplicity is different from theirs due to the presence of terms involving the NG functions 
El, £2, etc. in an exponential (see e.g. Eqn. (59)). Our prediction for the void multiplicity (or number density) 
is also simple to implement numerically, especially when combined with the fits for the NG functions given in 
Ref. [51]. 

It is worth spending a moment to compare the relative merits of using voids as opposed to halos as a probe of 
primordial NG. In terms of absolute numbers, a given survey (or simulation box) will always contain fewer voids 
of some Lagrangian scale R than halos of the corresponding mass scale M oc R 3 . This is simply because the 
void barrier height S v is larger than the halo barrier height S c , leading to a stronger cutoff ~ e~^ Vv for the void 
abundance. Furthermore on comparing the typical masses of the largest clusters observed (M ~ 10 15 h~ 1 M so \) 
with the typical comoving sizes of the largest voids (i? C om ~ 50ft. -1 Mpc), one finds that the largest Lagrangian 
length scales probed by both halos and voids are roughly the same (i? ~ 25ft. _1 Mpc). In this sense voids would 
be a poorer statistical probe of primordial NG than halos. The strength of voids however comes from the 
complementarity which is demonstrated in Fig. 2, and was also highlighted in Refs. [47, 48]. In contrast to 
the halo abundance which is e.g. reduced compared to the Gaussian for negative /nl, the void abundance is 
enhanced at large radii, and vice versa for positive /nl- R will be interesting to see how these characteristics 
ultimately play out in determining the constraining power of voids. 

Our work can be extended in more than one direction. Firstly, our arguments regarding the terms with 
large \j\ in Eqn. (68) were somewhat qualitative, although we expect them to be robust. In a future work, we 
will test these arguments by numerically generating the appropriate random walks and explicitly determining 
the multiplicity using the resulting distribution of conditional first crossing times [56]. Of course it will also 
be interesting to compare our predictions (and those of others) with full-fledged iV-body simulations. On the 
analytical side, it will be interesting to study the regime v y > 1 in more detail, since apart from void-in-cloud 
effects it is also likely that triaxial effects would become important here [57, 58]. It would be interesting to 
try and account for such effects (even approximately) within our framework, perhaps along the lines explored 
in Ref. [59] (see also Refs. [60, 48]). Finally, an equally interesting avenue would be to convert our predictions 
which hold for the dark matter distribution, to predictions for voids in the galaxy distribution, perhaps by 
generalizing the treatment of Ref. [46] to the non-Gaussian case. 
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Appendix 

A Primordial non-Gaussianity 

The physics of inflation governs the statistics of the initial seeds of inhomogeneities which grow into the 
large scale structure we see today. These initial conditions can be characterised by the primordial comoving 
curvature perturbation lZ(x) (with Fourier transform TZ(k)) which remains constant on superhorizon scales. 
The function {TZ(ki)TZ(k 2 )TZ(k 3 )) c (with the subscript denoting the connected part) is then an example of a 
function which probes the physics of the inflationary epoch. By translational invariance, it is proportional to 
a momentum-conserving delta function: 

(ft(ki)ft(k 2 )ft(k3) >c - (27r) 3 Mk! + k 2 + k 3 )B n (k u k 2 , h) , (A.l) 

where the (reduced) bispectrum B-R,(ki, fc 2 , fc 3 ) depends only on the magnitude of the k's by rotational invari- 
ance. According to the particular model of inflation, the bispectrum will be peaked about a particular shape 
of the triangle. The two most common cases are the squeezed (or local) NG, peaked on squeezed triangles 
fci <C fc 2 ~ fc 3 , and the equilateral NG, peaked on equilateral triangles ki ~ fc 2 ~ fc 3 . Indeed, one can define 
a scalar product of bispectra, which describes how sensitive one is to a NG of a given type if the analysis is 
performed using some template form for the bispectrum. As expected, the local and equilateral shapes are 
approximately orthogonal with respect to this scalar product [61]. We will now describe these two models in 
more detail. 

The local model: 

The local bispectrum is produced when the NG is generated outside the horizon, for instance in the curvaton 
model [62, 63] or in the inhomogeneous reheating scenario [64]. In these models, the curvature perturbation 
can be written in the following form, 

ft(x) = ft fl (x) + (^(x) -(H 2 g )) + ^Snl^(x) , (A.2) 

where lZ g is the linear, Gaussian field. We have included also a cubic term, which will generate the trispectrum 
at leading order. The bispectrum is given by 

Bn(h,k 2 , k 3 ) = If™ [P n {k 1 )Pn{k2) + cycl.] , (A.3) 

5 

where "cycl." denotes the 2 cyclic permutations of the wavenumbers, and Pn{k) is the power spectrum given 
by Pn{k) — Ak ns ~ 4 . The trispectrum is given by 



(7e(k 1 )^(k 2 )^(k 3 )^(k 4 ) ) c = (2 7 r) 3 fe(k 1 + k 2 + k 3 + k 4 ) 

• (A.4) 



36 54 

/^^P K (|k a +k 6 |)P R (y^(fc c ) + -3 NL Pn(ka)Pn{h)Pn(k c ) 



25 -^ ^ - »■/ "> w 25 

b<c a<b<c 
a^b,c 



The equilateral model: 

Models with derivative interactions of the inflaton field [65, 66, 67] give a bispectrum which is peaked around 
equilateral configurations, whose specific functional form is model dependent. Moreover, the form of the 
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bispectrum is usually not convenient to use in numerical analyses. This is why, when dealing with equilateral 
NG, it is convenient to use the following parametrization, given in Ref. [68], 
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Bate, fa, **) = y/N q L UI A + 3(fclfc2fc3)2 (4- n3 )/3 - (fc^^^/a + 5 P«™J ■ ( A - 5 ) 

This is peaked on equilateral configurations, and its scalar product with the bispectra produced by the realistic 
models cited above is very close to one. Therefore, being a sum of factorizable functions, it is the standard 
template used in data analyses. 

To connect the statistics of 1Z with large scale structure, we use the fact that the excursion set ansatz only 
requires us to know the linearly extrapolated present day behaviour of the density contrast 5r smoothed on 
various length scales. We can relate this quantity to the initial conditions 1Z via the following relations. We 
start from the Bardeen potential $ on subhorizon scales, given by 

9{k,z) = ~T(k)^-Tl(k), (A.6) 
o ct 

where T{k) is the transfer function of perturbations, normalized to unity as k — > 0, which describes the 
suppression of power for modes that entered the horizon before the matter-radiation equality (see e.g. Ref. [69]); 
and D(z) is the linear growth factor of density fluctuations, normalized such that D{z) = (1 + z)^ 1 in the 
matter dominated era. Then, the density contrast field is related to the potential by the cosmological Poisson 
equation, which in Fourier space reads 

Ink 2 Ik 2 

= M(k, z)TZ(k) , (A.7) 

with fl m the present time fractional density of matter (cold dark matter and baryons), and H = 100ft- km s~ Mpc~ 
the present time Hubble constant. The redshift dependence is trivially accounted for by the linear growth factor 
D(z). Introducing a filter function W.r(|x|), the smoothed density field (around one point, which we take as 
the origin) is given by 

d 3 k 
(2nf 



Sr(z)= I 7 — E W(kR)8(Kz), (A.8) 



where W(kR) is the Fourier transform of the filter function. The results of this paper are strictly valid only 
for a sharp filter in fc-space, although a physically more relevant filter would be the spherical top-hat filter in 
real space, whose Fourier transform W(kR) is given by 

3 

w(y) = ^ ( sin y - v cos y) ■ ( A - 9 ) 

While this choice allows us to have a well-defined relation between length scales and masses, namely M = 
(47r/3)fi m p c i? 3 with p c = 3H$/(8wG) = 2.77 • lO^^-Mso^/i^Mpc) -3 , it spoils the Markovianity of the 
random walk of Sr (see e.g. Ref. [52]). We will therefore continue to present results for the sharp-k filter. By 
using Eqns. (A.8) and (A.7) we have, for the 3-point function, 

/d 3 ki rl 3 ki d 3 ki ■ ■ 
(»3 (27T)3 (2 ^ ) 3 ^( fc l^)^( fc 2^2)^(fc 3 i 1 ' 3 )Al(fc 1 )Al(fc 2 )Al(fc3)(^(k 1 )^(k 2 )^(k 3 )) e , 

(A.10) 

where we suppressed the redshift dependence, and analogous formulae are valid for the higher order correlations. 
Fig. 5 shows the behaviour of the NG functions e\ and e 2 defined in the text (Eqn. (44)) as a function of S = a R 
in the local and equilateral cases. (e 2 is shown only for the local case.) We see that these functions remain 
approximately constant over a range of S values which corresponds to roughly the range 2-25/i _1 Mpc of 
Lagrangian smoothing length scales. 
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Figure 5: Scale dependence of the e n . Panel (a) : Behaviour of e± vs. a 2 in the local and equilateral models, for 
/nl = 100 in each case. Panel (b) : Behaviour of e-i for the local model with /nl = 100 and <7nl = 10 4 . The terms 
proportional to / N l an d <?nl are shown separately. Also shown is ef for the same model. The axes are logscale. 

B Some results concerning statistics of random walks 
B.l Expressing J^vdW as a sum of Gaussians 

Let J- denote the r.h.s of Eqn. (8). To prove the equivalence of the two expressions for J-"svdw i n Eqns. (5) and 
(8), we will show that the Laplace transform C{uj) of T is the SvdW result (6). We have 

C{u) = / &Se-" s F{S) 
Jo 



E 



S v - 2jS T 



+ 00 



dS 



(2tt) 1 /2 J q S 3 ^ 
]T sgn(2^ T - Sy) e -IW*r-«rl^ 



-uiS-(5 v -2j8 T ) 2 /2S 



j=—oo 



j<0 j>0 



(B.l) 



w 



here the integral in the second line can be calculated using Eqn. 3.472(5) of Ref. [70], and we used the fact 
that 2j<5t — 5 V is negative for j < 0. By changing j — > —j in the first summation in the last line one gets 



C(cj)=e- S ^Y1 

3=0 



,-25tV2w 



+00 

; E 

7 = 1 



,-2<5tV5w 



-S v \/2uj _ g(5 v V2"g-2i5TV / 2aJ 
^ g— 2S t V2lJ 



g<5 c g— 5 c v / 2uJ 

g^T y/2oj g— St ^/2ui 



(B.2) 



which reproduces exactly the result of Eqn. (6). 



B.2 Wrong side counting for the single barrier 

Here we show that the wrong side (WS) counting argument for the single barrier, is equivalent to MR's 
calculation of the fx. rate for generic random walks. For a single barrier problem with discrete time steps, the 
probability for the walk to remain below the barrier S c for the first n — 1 steps, and to cross the barrier at the 



,th 



step is dS n J 5 ^ dS n -i . . . ddi W(6 ; {Sk} n ; S n )- Here W(5 ; {5k} n ; S n ) defined in Eqn. (14) is the p.d.f. 
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of a completely generic random walk - we are not assuming that the walk is Markovian or even Gaussian. 
Using an argument very similar to the one used in writing Eqn. (23), the WS f.c. rate is then given by 

-7 r ws,i-bar(<S') = lim —— / dS n d6 n -i . . .dS 1 W(S ;{S k } n ; S n = S) . (B.3) 
as^o AS J 5c J_ oa 

In general we also have the identity 

/oo pS c p8 c 

dS n d5 n - 1 ...d5 1 W(5 ;{6 k } n] S n )= dS^ . . . dW(<5 ; {6 k } n -i; S n -i) , (B.4) 
-oo J —oo J —oo 

which is a simple marginalisation over the location of the walk at the final time step, and follows from the 
definition of W(8q; {Sk} n ', S n ). Let Pas(> S n ) denote the probability that the f.c. time is larger than S n . 
Since this is the same as the probability that the walk has not crossed the barrier in the first n steps, we have 

Pas(> S n ) = ( C dS n ... dd^iSo; {5 k } n ; S n ) , (B.5) 



and the r.h.s. of Eqn. (B.4) equals Pas(> S n -\). Splitting the integral over S n on the l.h.s. of Eqn. (B.4) as 
— » f ^ + and using Eqn. (B.5) we then get 

/'OO rS c 

/ dS n dS n - 1 ...d5 1 W(5o;{5k} n ;S n ) = P A s(>S n - 1 )-P AS (>S n ). (B.6) 

J 5 C J — oo 

From Eqn. (B.3) it then follows that 

^WS,l-bar(S) = lim o [P AS (> S„_i) - P AS (> S n )} = -d S PAS=0(> S) = J- M R,l-bar(S) , (B.7) 

where Pas=o(> S) is the continuum limit cumulative probability for the f.c. time, which is what MR use to 
compute the f.c. rate. This proves our result. 

B.3 Comparison with Lam et al. [48] 

To compare with the analysis of Lam et al. [48] it is useful to express the two barrier conditional f.c. rate in a 
slightly different form. First, by manipulating the integrals in its definition, the conditional probability density 
nAs(<5o, 5 n ; S) can be rewritten as 

^■As{So,S n ;S n ) = J dS n -i . . .dSiW(d ;{Sk} n ;S n ) 
dSn-i ■ ..dSiW(So; S n ) 

rl ~_ roo />oo />5 C 

- ^ / d5 n -i . . . d5 k+ i / dd k / dS k -i...dSiW(S ;{S k } n ;S n ), (B.8) 

t — -i J — 5 V J S c J — 5 V 



/oo 



k =l ' 



where it is understood that for k = 1 the integral in the second term begins with d8\ on the far right, and 
for k — n — 1 it ends with dS n -i on the far left. This form of the expression is useful because it helps us 
to intuitively understand the calculation of the conditional f.c. rate, which becomes 



Fng{S) = lim — / d5 n U A s(So,S n ;S) 
as^o AS J_ oa 

= lim ITS / d<5 « / d5 n - 1 ...d5 1 W(5 ;{5 k } n ;S) 

AS^O AS J_ oa J_x 



1-8 

n-1 r S 



Y /"-Ov roo f>oo />o c 

Jjp 1 l^d^2 dS n dSn-! . . . d5 k+1 / dS k / dS k _ 1 ...dS 1 W(S Q ;{d k } n ;S). 
as^o AS J_ x J_ Sv J Sc J_ 5v 

(B.9) 
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The first term here is simply the single barrier rate across the void barrier (see Appendix B.2). The second term 
subtracts from this rate the fraction of trajectories which crossed the halo barrier at times before S, without 
having crossed the void barrier before. This is exactly what Lam et al. write in Eqn. 41 of Rcf. [48] (in the 
continuum limit the summation over k above would become an integral over time) . Their expression however 
assumes a factorisation of the second term into a two barrier rate convolved with a conditional probability. 
This is certainly true for the Gaussian case with the sharp-k filter, as one can immediately see by replacing W 
above with W gra which factorises, giving the two barrier probability across the halo barrier, convolved with the 
usual Gaussian conditional probability. In the non-Gaussian case however, W(Sq; {Sk} n ] S) does not factorise 
due to the presence of multi-scale correlations which make the process non-Markovian. In practice though 
these effects are small, and hence factorisability should be a reasonable assumption. Notice however that our 
approach does not require this assumption, and we also do not need to linearize the effects of non-Gaussianities 
as in Ref. [48]. 

C Details of the calculation of J-"s v dW 

C.l T^ s {5) for fixed 6 < -S v and small AS 

With S < — S v held fixed, a few manipulations allow us to write I^ s (5) as 



rAS 



(6 < -5v)) 



(2A5) 



i/2 



7rV2 

(2A5)"/ 2 



7T 



1/2 



-iy 



(<5+<5 v )/v / 2A5 



{S-S C )/V2AS 



dyy n e- y 



dyy n e- y 



\S+5 V \/V2AS 



(C.l) 



where we ignore the second term since one can explicitly check that it leads to an exponentially suppressed 
contribution which eventually vanishes in the appropriate continuum limit. (This essentially follows from the 
fact that St = S c + S v is finite and hence 8t/V2AS — > +oo as AS — > 0.) Now, since the variable q in the last 
integral is strictly positive under our assumptions, one can check that the transformation to x = q 2 will reduce 
the integral above to an incomplete Gamma function, 



Z£ S (5 < Sy) 



(2AS) 



i/2 



1-1/2 



-(-!)' 



ir 

2 



n + 1 (S + 5 V ) 2 



2AS 



(C.2) 



and in the limit of small AS, the asymptotic expansion of the incomplete Gamma gives us Eqn. (29). 

D Details of the calculation of Jng 

D.l The two barrier derivative exchange property 

To prove Eqn. (51), for any function g{8\, . . . , S n ) define the function h(5 n ;S c ,S v ) of three variables as the 
multiple integral 



h(S n ; S c , S v ) = J dSn-i ...d5\ g{5\, ...,$„). 
The l.h.s. of Eqn. (51) can then be reduced as 



(D.l) 



/ — d v rd c _ _ 

dS n / dS n -i . . . dS\ djg 
-oo J — 5 V 

n-l p s v P S C 



= h(-5 v ]5 c ,S Y ) + / dS n I dS n -i...dSidjg 

J —oo J — 5v 

= h(-5 v ;5 c ,S v ) + / dd n (d Sc \s v - ds v \s c ) h{S n \ S C ,S V ) . 

J —oo 



(D.2) 
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where in the second line we integrated the term involving d$ n and the third line follows since for each 1 < j < 
n — 1, the integral becomes 



/— <5 V ffi c 
d5„ I d5 n -i . . . 
-co J — S v 



dSj + idSj-i . . 



.dS 1 



l(8i, ...,6j = 6 C ,.. .,8 n ) - g(6i, ...,Sj = -S v , ...,8 n ) 



(D.3) 



the summation of which is the same as the integral in the third line. The derivative d$ c \s v now simply comes 
across the integral over 8 n since the boundary is independent of 8 C . The term involving d$ v \s c can be handled 
by noting that when this derivative acts only on the boundary of the integral J_ v dS n h(5 n ; 8 C , 8 V ), we simply 
get — h(— 8 V ; 8 C , 8 V ). We can then write 



d sJs. 



— <5 V p — <5 V 

d8 n h(8 n ;8 c ,S v ) = -h(-8 v ;8 c ,8 v ) + / dS n ds v \s c h(S n ;S c ,S v ) . 
-00 J —00 



This simplifies the expression (D.2) to 



L.h.s. = (<9«5 C |,5 V - d Sv \s, 



J — < 



d5 n h(5 n ; S Cl S v ) . 



(DA) 



(D.5) 



It is also straightforward to show that under the change of variables (S c , S v ) — > (St, S v ) where St = S c + 8 V , we 
have ds c \s v - ds v \s c = -dsJs T , which gives us 



L.h.s. = - 



d 



5t 



/—S v p8 c 
dS n / 
-00 J — s v 



d5 n -i ...dSi g(5i, . . . , 5 n ) = R.h.s. , 



which proves the result. 



(D.6) 



D.2 Leading unequal time contribution 

Here we sketch a proof of Eqn. (54) . Defining the quantity Mj as 

-<5 V rSc 

we can write 

j 

Using the derivative exchange property (51) we can write Mj as 



. 1_1 J— 00 J— I 



fe,/=l 



m, = at 



dS n 



d5 n ^ 1 ...d8 1 d J W &m 



(D.7) 
(D.8) 



dS n pS(<5o, ~8 V ; Sj)U^ s (-8 v , 8 n ; S - Sj) - U s ™ s (8 , S c ; ^)n| m s (<5 c , S n ; S - Sj)] . (D.9) 



Recall that ds v = d/dS v \s T . In the second line, we already have an expression for 11^ (do, — S v ; Sj) from 
Eqn. (38). One can also check that in the continuum limit, the second term involving S c will be exponentially 
suppressed. The only non-trivial quantity we need then is the integral J_ Sv dS„H^(—S v , 5 n ;S — Sj). Following 
techniques similar to those discussed in Section 3.2 leading up to Eqn. (34), we can write this integral by simply 
replacing 5q in (34) with (— S v ) to get 



dd n U s ^ a s (-d v ,d n ; S - S^) = 



- ( AS 



V 2tt 



1/2 



drj 



v(rj)e n 



n^ s .(— 5 V , — S v ; S Sj) , 



(D.10) 



with the same boundary layer function v(rj) appearing, since this is independent of the value of 5q. Notice 
that the integral involving v(rj) is identical to that appearing in the expression (39) for Jsvdw- In fact the 
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other constant 7 in that expression also appears in Mj through the object (So, —S v ; Sj), and will lead us 
to identify a factor of Jsvdw in the final expression, see below. The object n^(— S v , — S v ; S — Sj) is similar to 
one discussed by MR in Ref. [52], where they showed that n^ rn bar (S C ,S C ;S) = AS/ (2ns 3 ) 1 / 2 , independent 
of <5 C . In fact, one can show that H^(— <5 V , — S V :S) is also independent of both S c and S v , and has the same 
value 



U^ S (-S V -S V ;S) 



AS 1 



(D.ll) 



/2nS 3 / 2 

This can be checked in two steps. Firstly for arbitrary n one can show that the path integral with n steps 
involved in n^(— S v , — <5 V ; S n ) becomes independent of S c and S v in the continuum limit (since S c and S v only 
appear in the limits of integration, and a suitable change of variables sends these to infinity as AS — > 0). 
The remaining integral is therefore only a function of n. Following MR, a dimensional argument coupled with 
a straightforward calculation for n = 2 then fixes the time dependence completely, leading to Eqn. (D.ll). 
Putting things together, 



Mj = -dl 



n^ s (S ,-S v ;S j )i — 




drj 

i-v) 



v(rj)e 




1 



^(S-Sjf/ 2 



(AS) 



2 FsvdW (Sj) 

2 (S-S,) 3 / 2 



■dl 



(D.12) 



where we used Eqns. (38) and Eqn. (39) in writing the last line. Using this result for Mj in Eqn. (D.8) with the 
continuum limit £"=1 h(Sj) -> J dSh(S)(AS)- 1 , we see that the factors of AS cancel, and further expressing 
Q^'°'°\S) in terms of c±(S) finally leads to the result of Eqn. (54). 
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